A microbial diversity perspective within the Restoring Resilient Ecosystems project

Dig Deeper - Microbial Community Analysis

Introduction

Biodiversity loss, ecosystem degradation, and habitat destruction are increasingly linked to human-driven changes in land use, including urbanisation, agriculture, and the exploitation of natural resources (European Parliament, 2025; Jaureguiberry et al., 2022). In response, governments across Europe — including the EU — have introduced ambitious environmental strategies such as the EU Biodiversity Strategy for 2030 (European Parliament, 2025) and the 30x30 target (Markwick, 2023), which aims to protect 30% of land and sea by the year 2030.

Ecological restoration plays a vital role in addressing these challenges. Rather than simply returning ecosystems to a previous state, modern approaches focus on restoring ecological processes and enhancing ecosystem resilience (Hicks, 2023).

The RestREco Initiative

RestREco (Restoring Resilient Ecosystems) is a NERC-funded research project that adopts a resilience-based perspective on ecological restoration. The initiative brings together researchers from:

• Cranfield University
  
• University of Stirling
• UK Centre for Ecology & Hydrology
  
• The National Trust
  
• Forest Research

Using a natural experiment design, RestREco studies a network of 133 ecological restoration sites across England and Scotland. The project aims to identify key drivers of ecosystem development, such as:

• Time since restoration began
  
• Initial ecological conditions
  
• Proximity to existing woodland and grassland

The goal is to understand how these factors influence ecosystem complexity, function, and resilience to future pressures (RestREco, 2024).

The Dig Deeper Study

As part of the RestREco initiative, the Dig Deeper study focused on how the age of restoration, establishment type, and site management affect soil microbial communities, specifically bacteria and fungi.

To explore this, high-throughput sequencing was conducted on:

• 16S rRNA gene (for bacterial communities)
  
• ITS region (for fungal communities)

The analysis focused on three main aspects:

• Alpha and beta diversity
  
• Taxonomic composition
  
• Functional diversity

These microbial assessments complement broader ecosystem-level measurements within the RestREco project, including vegetation, invertebrates, and ecosystem functions such as litter decomposition, pollination services, and soil thermodynamic efficiency.

The following sections describe the sampling design, metadata structure, and the processing pipeline used to characterise microbial communities.

---

Sample Design and Metadata Overview

Sample Collection and Geographic Coverage

A total of 330 soil samples were collected in 66 sites of England for each marker (5 per site).

Sample Zone - Based on GPS coordinates

Overview of the soil sampling and sequencing strategy for each microbial marker.

Metric	16S	ITS
Microbial group	Bacteria	Fungi
Region sampled	England	England
Number of sites	66	66
Samples per site	5	5
Total samples	330	330
Average reads per sample	~65,000	~65,000
Read count range	30,000–85,000	10,000–90,000

Sampling Summary

Metadata Overview

Each sample collected was accompanied by metadata capturing key environmental and management variables. These contextual factors were essential for interpreting variation in microbial diversity.

Description of metadata variables associated with soil samples

Variable	Description
Site	Name of the sampling site
Plot number	Subdivision of each site (usually 5 plots per site)
CU Code	Unique code for each sample
Year_est	Year of establishment of the site
Age	Site age (ranging from 1 year to over 100 years )
Latitude/Longitude	GPS coordinates of the sample
Establishment	Restoration type or land management
pH	Soil pH value at the time of sampling
EC	Electrical conductivity of the soil
Cutting	Whether the site is cut (1 = Yes, 0 = No)
Cattle	Presence of cattle grazing (1 = Yes, 0 = No)
Sheep	Presence of sheep grazing (1 = Yes, 0 = No)
Plough	Whether the soil has been ploughed (1 = Yes, 0 = No)

Metadata Summary

---

Analysis Pipeline

Input Files Selection

Before proceeding with bioinformatic analyses, it was necessary to select consistent input files for each sample. For both ITS and 16S datasets, each sample was associated with multiple types of sequencing files. For example, a sample such as GF677 had up to five different files: GF677_1.fastq.gz, GF677_2.fastq.gz (paired-end reads with barcode and primer removed), GF677.raw_1.fastq.gz, GF677.raw_2.fastq.gz (raw unprocessed reads), and GF677.extendedFrags.fastq.gz (merged forward and reverse reads). To ensure consistency and avoid redundancy, we selected one file type per sample for downstream analysis. For the 16S dataset, we used the paired-end reads with barcode and primer removed (*_1.fastq.gz and *_2.fastq.gz), while for the ITS dataset, we selected the raw reads (*.raw_1.fastq.gz and *.raw_2.fastq.gz), which were best suited for our quality filtering and denoising steps.

Denoising and Feature Table Construction

The sequencing data were processed using the QIIME 2 bioinformatics platform — a widely used tool for microbiome analysis. Raw amplicon reads were denoised using the DADA2 plugin, enabling accurate identification of amplicon sequence variants (ASVs) with single-nucleotide resolution. This step also included chimera removal and read quality trimming, ensuring high-confidence input for downstream analyses.

Following denoising, a feature table was constructed for each dataset (16S and ITS), summarising the number of sequences associated with each ASV across samples.

Diversity Analysis and Taxonomic Assignment

Following quality control and feature extraction, downstream analyses were performed to characterise microbial diversity and taxonomic composition.

Alpha and beta diversity metrics were computed using the diversity and emperor plugins in QIIME 2, enabling comparison of microbial communities across restoration gradients and site conditions.

Taxonomic classification of ASVs was then performed using trained classifiers against reference databases: GreenGene2 for bacterial 16S sequences and UNITE for fungal ITS sequences. This allowed each ASV to be annotated with its likely taxonomic lineage (from kingdom down to genus or species when possible).

TO BE COMPLETED/CHANGED

Summary of Pipeline

Workflow (16S)

---

Data Pre-Processing

MutliQC on raw data

Bacteria (16S)

You can explore the full MultiQC report by clicking the image below:

MultiQC Plot (16S)

Fungi (ITS)

You can explore the full MultiQC report by clicking the image below:

MultiQC Plot (ITS)

QC after denoising

Statistics Table Summary

After denoising, quality control (QC) metrics provide essential insights into the effectiveness of the data processing steps and the overall quality of the resulting feature table. This summary table presents key statistics for each sample, including the number of input reads, filtered reads, and final feature counts after denoising. These metrics help assess sequencing success, identify potential outliers, and ensure that sufficient data remain for robust downstream analyses. Samples with unusually low read counts or feature richness may need to be excluded or interpreted with caution.

Bacteria (16S)

Here is a link to the statistics after denoising to view it on QIIME2 (16S) : Statitics after denoising (16S)

Statitics after denoising (16S)

Fungi (ITS)

Here is a link to the statistics after denoising to view it on QIIME2 (ITS) : Statitics after denoising (16S)

Statitics after denoising (ITS)

QC Plots

These Quality Control (QC) plots were generated after trimming the sequencing reads. They provide a visual summary of the base quality scores, read length distributions, and other metrics, helping to assess whether the trimming step successfully removed low-quality regions and adapter contamination. Consistently high-quality reads across samples are essential for reliable downstream analysis.

Bacteria (16S)

QC plot (16S)

Fungi (ITS)

QC plot (ITS)

Rarefication Curves

Rarefaction curves provide a visual tool to assess sequencing depth and compare species richness between samples. In this study, rarefaction curves were generated separately for bacterial (16S rRNA gene) and fungal (ITS) communities, using the number of observed features—i.e., unique ASVs—as a proxy for richness.

For both bacteria and fungi, the shape of each curve indicates whether sequencing depth was sufficient to capture most of the diversity in a given sample. Curves that level off suggest that a representative portion of the community has been sampled, whereas rising curves indicate that additional sequencing could reveal further diversity.

This step is crucial to ensure that downstream diversity analyses are not biased by unequal sampling effort.

Bacteria (16S)

Here is a link to the Rarefiction plots for more flexibility on QIIME2: Rarefiction plots (16S)

Select variable: Establishment Age Year of Establishment Os Location Plough Lat/long

Rarefaction Curves of Observed Features by 'Establishment' (16S)

Rarefaction Curves of Observed Features by 'Age' (16S)

Rarefaction Curves of Observed Features by 'Year_est' (16S)

Rarefaction Curves of Observed Features by 'OS_location' (16S)

Rarefaction Curves of Observed Features by 'Lat_long' (16S)

Rarefaction Curves of Observed Features by 'Plough' (16S)

---

Fungi (ITS)

Here is a link to the Rarefiction plots for more flexibility on QIIME2: Rarefiction plots (ITS)

Select variable: Os Location CU Code Lat/long

Rarefaction Curves of Observed Features by 'OS_location' (ITS)

Rarefaction Curves of Observed Features by 'CU Code' (ITS)

Rarefaction Curves of Observed Features by 'Lat_long' (ITS)

---

Taxonomic Diversity (Bacteria & Fungi)

Alpha Diversity

Alpha diversity refers to the variety of organisms within a particular sample or environment. It reflects both richness—the number of distinct taxa—and evenness—how evenly individuals are distributed among those taxa. One of the most widely used measures for assessing alpha diversity is the Shannon index.

The Shannon index takes into account not only the number of species present, but also how evenly their abundances are distributed. A higher Shannon value generally indicates a more diverse and ecologically balanced community.

Another important metric is Faith’s Phylogenetic Diversity (Faith PD), which measures the total branch length of the phylogenetic tree that spans the species in a sample. Unlike the Shannon index, Faith PD incorporates evolutionary relationships, providing a phylogenetic perspective on diversity.

We also include Pielou’s Evenness index, which specifically quantifies how equally individual organisms are distributed across taxa. While Shannon integrates both richness and evenness, this metric isolates the evenness component, providing a complementary view of diversity patterns.

To allow interactive exploration of alpha diversity metrics across different environmental variables, we implemented a drop-down menu that dynamically displays the corresponding plots. Some variables, such as pH category, are only present in the ITS dataset, while others, like Year group, are specific to the 16S dataset. Internally, variables are mapped to their dataset-specific equivalents where needed (e.g. Age group in 16S becomes Age category in ITS). It is important to note, however, that these variables are not always directly comparable: for instance, Age group (16S) divides sites into multiple discrete intervals based on restoration age, while Age category (ITS) is a binary classification based on whether a site is above or below the median age. Despite these differences, the interface ensures that only available and relevant plots are shown for each selection.

Bacteria (16S)

In the plots below, we examine how the Shannon index, Faith PD and Evenness vary across different environmental and experimental conditions for the 16S (bacteria).

Select variable: Establishment Cutting Cattle Sheep Plough Year group Age group pH binary Age binary

Shannon Index Boxplot

The boxplot below illustrate differences in Shannon diversity across groups. This metric reflects both species richness and how balanced the community is in terms of species abundance.

Here is a link to the full QIIME2 results (16S) : Shannon Index (16S)

Shannon entropy revealed significant overall differences, with green hay/bush sites showing higher diversity than both natural regeneration and seed mix sites. Seed mix sites also had higher Shannon diversity than natural regeneration.

Kruskal-Wallis p-value: 0.00102

Kruskal-Wallis p-value: 0.158

Kruskal-Wallis p-value: 0.939

Kruskal-Wallis p-value: 0.655

Kruskal-Wallis p-value: 0.316

Kruskal-Wallis p-value: 0.408

Kruskal-Wallis p-value: 0.126

Kruskal-Wallis p-value: 0.116

Kruskal-Wallis p-value: 0.358

Faith PD Boxplot

The following plots show Faith’s Phylogenetic Diversity, which integrates evolutionary relationships to capture how phylogenetically broad each microbial community is.

Here is a link to the full QIIME2 results (16S) : Faith PD (16S)

Pielou’s evenness differed significantly across establishment methods (p = 0.0148), with green hay/bush sites being significantly more even than natural regeneration (q = 0.0148). The difference between seed mix and natural regeneration approached significance (q = 0.0597), while no significant difference was found between green hay/bush and seed mix.

Kruskal-Wallis p-value: 0.0106

Kruskal-Wallis p-value: 0.362

Kruskal-Wallis p-value: 0.878

Kruskal-Wallis p-value: 0.604

Kruskal-Wallis p-value: 0.705

Kruskal-Wallis p-value: 0.216

Kruskal-Wallis p-value: 0.0709

Kruskal-Wallis p-value: 0.0893

Kruskal-Wallis p-value: 0.393

Evenness Boxplot

These boxplots display Pielou’s Evenness, highlighting how uniformly taxa are represented in each community. It allows us to isolate imbalance in dominance from richness effects.

Here is a link to the full QIIME2 results (16S) : Pielou Evenness (16S)

Faith’s PD also varied significantly by establishment method (p = 0.0106), with both green hay/bush (q = 0.0148) and seed mix sites (q = 0.043) exhibiting greater phylogenetic diversity than natural regeneration. No significant difference was observed between green hay/bush and seed mix.

Kruskal-Wallis p-value: 0.0148

Kruskal-Wallis p-value: 0.17

Kruskal-Wallis p-value: 0.878

Kruskal-Wallis p-value: 0.291

Kruskal-Wallis p-value: 0.64

Kruskal-Wallis p-value: 0.854

Kruskal-Wallis p-value: 0.608

Kruskal-Wallis p-value: 0.542

Kruskal-Wallis p-value: 0.0817

---

Fungi (ITS)

For the ITS dataset, the alpha diversity was done per sample. In order to align it with the 16S analysis—where samples were already grouped by site—we aggregated the alpha diversity values by computing the mean per site. Categorical metadata was simplified using the most common (modal) value per site. This ensures consistency across datasets in the visual outputs. However, users interested in the original, unaggregated sample-level data can explore the full QIIME 2 results via the links provided under each section.

For the variable pH_category, the median is 7.94 and for Age_category it’s 14 years.

Select variable: Establishment Cutting Cattle Sheep Plough Age category pH category

Shannon Index Boxplot

The boxplot below illustrate differences in Shannon diversity across groups. This metric reflects both species richness and how balanced the community is in terms of species abundance.

Here is a link to the full QIIME2 results (ITS) : Shannon Index (ITS)

Fungal alpha diversity did not differ significantly between establishment methods or management types. However, Shannon diversity was significantly lower in alkaline soils (pH > 7.9, p < 0.02) and significantly higher in older sites (age > 14 years, p < 0.02), indicating that both soil pH and site age influence fungal richness and evenness.

Kruskal-Wallis p-value: 0.796

Kruskal-Wallis p-value: 0.719

Kruskal-Wallis p-value: 0.768

Kruskal-Wallis p-value: 0.308

Kruskal-Wallis p-value: 0.647

Kruskal-Wallis p-value: 0.159

Kruskal-Wallis p-value: 0.193

Evenness Boxplot

These boxplots display Pielou’s Evenness, highlighting how uniformly taxa are represented in each community. It allows us to isolate imbalance in dominance from richness effects.

Here is a link to the full QIIME2 results (ITS) : Pielou Evenness (ITS)

Sheep grazing was the only factor significantly affecting fungal evenness (p = 0.03), with grazed sites showing slightly higher Pielou’s Evenness. No significant effects were observed for establishment method, site age, or soil pH.

Kruskal-Wallis p-value: 0.512

Kruskal-Wallis p-value: 0.719

Kruskal-Wallis p-value: 0.928

Kruskal-Wallis p-value: 0.0927

Kruskal-Wallis p-value: 0.389

Kruskal-Wallis p-value: 0.386

Kruskal-Wallis p-value: 0.954

---

Comparative Microbial Community Composition (Beta Diversity)

To explore differences in microbial communities, we often rely on dimensionality reduction techniques such as Principal Coordinates Analysis (PCoA), visualised through Emperor plots. Two commonly used distance metrics in this context are Bray-Curtis and Jaccard.

While both metrics can reveal meaningful clustering and separation in microbial data, they capture complementary aspects of community structure.

Bray-Curtis

Emperor Plot

The Bray-Curtis Emperor plot is a 3D visualisation of microbial community dissimilarities between samples, based on the Bray-Curtis distance. This distance metric quantifies how different two samples are in terms of species abundance, taking into account both presence/absence and relative abundances. It does not incorporate evolutionary relationships between features.

Using Principal Coordinates Analysis (PCoA), the high-dimensional Bray-Curtis distance matrix is projected into a lower-dimensional space—typically three axes—to capture the main patterns of variation across samples.

The Emperor plot is an interactive 3D tool developed for QIIME 2 that allows users to explore these PCoA results. Samples are represented as points in space, and their spatial proximity reflects ecological similarity:

• Samples that are closer together have more similar microbial communities.
• Samples that are further apart differ more strongly in community composition.

This type of plot is particularly useful for identifying clustering by experimental or environmental factors—such as establishment type, land management practices, or age group—and for detecting gradients or patterns in microbial community composition.

Bacteria (16S)

Bray-Curtis beta diversity analysis identified establishment method, pH, and sheep grazing as the strongest drivers of differences in bacterial community composition across sites. Site age also showed a significant effect with this metric.

PCoA plots based on Bray-Curtis distances revealed clustering of samples according to establishment method, with visible separation between natural regeneration, green hay/bush, and seed mix groups, although some overlap was observed.

Here is a link to the Bray-Curtis Emperor Plot for more flexibility on QIIME2: Bray-Curtis Emperor Plot (16S)

Bray-Curtis Emperor Plot

Fungi (ITS)

Bray-Curtis analysis identified establishment method, soil pH, and site age as the primary drivers of fungal community differences, followed by sheep grazing, ploughing, and cutting. In PCoA plots, natural regeneration sites formed a loose cluster near seed mix sites, while green hay/bush sites clustered more distinctly along Axis 1.

Here is a link to the Bray-Curtis Emperor Plot for more flexibility on QIIME2: Bray-Curtis Emperor Plot (ITS)

Bray-Curtis Emperor Plot (ITS)

Permanova using Bray-Curtis distance – Taxonomic Level

The Bray–Curtis distance takes into account the relative abundance of taxonomic features, making it suitable for detecting quantitative differences in microbial communities across groups. When used in a PERMANOVA, this metric tests whether the overall structure and abundance of taxa vary significantly according to explanatory variables such as establishment type, age of the grassland, year of establishment, or management practices including cutting, ploughing, and the presence of grazing animals.. This approach highlights shifts in dominant or highly represented taxa.

Bacteria (16S)

Bray-Curtis-based PERMANOVA revealed significant differences in bacterial community composition across establishment methods, with all pairwise comparisons between methods also significant. These results indicate that restoration strategy has a strong influence on both the composition and abundance structure of bacterial communities.

Site age also had a significant but weaker effect. However, clustering by age in Bray-Curtis PCoA plots was less pronounced, suggesting more subtle shifts in community structure over time.

Soil pH significantly shaped bacterial community composition, with Bray-Curtis PCoA plots showing clear clustering between high and low-to-neutral pH soils

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough Year group pH binary Age binary

View full QIIME2 results (Bray–Curtis – Establishment)

Figure: PERMANOVA Bray–Curtis for “Establishment” (16S)

View full QIIME2 results (Bray–Curtis – Cutting)

Figure: PERMANOVA Bray–Curtis for “Cutting” (16S)

View full QIIME2 results (Bray–Curtis – Cattle)

Figure: PERMANOVA Bray–Curtis for “Cattle” (16S)

View full QIIME2 results (Bray–Curtis – Sheep)

Figure: PERMANOVA Bray–Curtis for “Sheep” (16S)

View full QIIME2 results (Bray–Curtis – Plough)

Figure: PERMANOVA Bray–Curtis for “Plough” (16S)

View full QIIME2 results (Bray–Curtis – Year_group)

Figure: PERMANOVA Bray–Curtis for “Year_group” (16S)

View full QIIME2 results (Bray–Curtis – pH_binary)

Figure: PERMANOVA Bray–Curtis for “pH_binary” (16S)

View full QIIME2 results (Bray–Curtis – Age_binary)

Figure: PERMANOVA Bray–Curtis for “Age_binary” (16S)

Fungi (ITS)

PERMANOVA analysis identified pH category (p = 0.001) and establishment method (p = 0.002) as the strongest drivers of differences in fungal communities. Site age (p = 0.024) and sheep grazing (p = 0.025) also had significant effects.

Pairwise comparisons showed the strongest differences between green hay/bush and seed mix sites (p = 0.008), followed by green hay/bush vs. natural regeneration (p = 0.011). Natural regeneration sites exhibited the highest variation in community composition, with Bray-Curtis distances ranging from 0.25 to 0.95.

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough pH category Age category

View full QIIME2 results (Bray–Curtis – Establishment)

Figure: PERMANOVA Bray–Curtis for “Establishment” (ITS)

View full QIIME2 results (Bray–Curtis – Cutting)

Figure: PERMANOVA Bray–Curtis for “Cutting” (ITS)

View full QIIME2 results (Bray–Curtis – Cattle)

Figure: PERMANOVA Bray–Curtis for “Cattle” (ITS)

View full QIIME2 results (Bray–Curtis – Sheep)

Figure: PERMANOVA Bray–Curtis for “Sheep” (ITS)

View full QIIME2 results (Bray–Curtis – Plough)

Figure: PERMANOVA Bray–Curtis for “Plough” (ITS)

View full QIIME2 results (Bray–Curtis – pH_category)

Figure: PERMANOVA Bray–Curtis for “pH_category” (ITS)

View full QIIME2 results (Bray–Curtis – Age_category)

Figure: PERMANOVA Bray–Curtis for “Age_category” (ITS)

Jaccard

Emperor Plot

The Jaccard Emperor plot provides a 3D visualisation of microbial community dissimilarities based on the Jaccard distance. Unlike Bray-Curtis, the Jaccard metric considers only the presence or absence of features (e.g., microbial taxa) in each sample, ignoring their relative abundances.

This makes the Jaccard distance particularly suited for assessing community membership rather than abundance structure—focusing on which species are present, regardless of how abundant they are.

Using Principal Coordinates Analysis (PCoA), the high-dimensional Jaccard distance matrix is projected into a lower-dimensional space—usually three principal axes—to reveal major patterns in sample composition.

As with Bray-Curtis, the Emperor plot allows for interactive exploration of these ordinations:

• Samples positioned closely together share more taxa in common (i.e., similar membership).
• Samples far apart have fewer shared taxa, reflecting greater differences in species presence.

The Jaccard plot is useful when exploring factors that influence community membership, such as habitat type, land use, or environmental filtering—especially in studies where presence/absence patterns are more meaningful than relative abundances.

Bacteria (16S)

Here is a link to the Jaccard Emperor Plot for more flexibility on QIIME2: Jaccard Emperor Plot (16S)

Jaccard Emperor Plot (16S)

Fungi (ITS)

Jaccard analysis confirmed similar drivers of community structure—establishment method, pH, and age—with natural regeneration, seed mix, and green hay/bush sites forming distinguishable but overlapping clusters in higher dimensions. This suggests that both presence/absence and abundance-based differences in fungal communities are shaped by restoration strategy and site conditions.

Here is a link to the Jaccard Emperor Plot for more flexibility on QIIME2: Jaccard Emperor Plot (ITS)

Jaccard Emperor Plot (ITS)

Permanova using Jaccard distance – Taxonomic Level

The Jaccard distance considers only the presence or absence of features, offering a more qualitative perspective on microbial community differences. In this PERMANOVA, the focus is on whether the composition of taxa — regardless of their abundance — differs significantly between groups. This metric is particularly useful for identifying patterns in taxon occurrence, including rare or transient species that may not strongly influence abundance-based distances.

Bacteria (16S)

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough Year_group pH binary Age binary

View full QIIME2 results (Jaccard – Establishment)

Figure: PERMANOVA Jaccard for “Establishment” (16S)

View full QIIME2 results (Jaccard – Cutting)

Figure: PERMANOVA Jaccard for “Cutting” (16S)

View full QIIME2 results (Jaccard – Cattle)

Figure: PERMANOVA Jaccard for “Cattle” (16S)

View full QIIME2 results (Jaccard – Sheep)

Figure: PERMANOVA Jaccard for “Sheep” (16S)

View full QIIME2 results (Jaccard – Plough)

Figure: PERMANOVA Jaccard for “Plough” (16S)

View full QIIME2 results (Jaccard – Year_group)

Figure: PERMANOVA Jaccard for “Year_group” (16S)

View full QIIME2 results (Jaccard – pH_binary)

Figure: PERMANOVA Jaccard for “pH_binary” (16S)

View full QIIME2 results (Jaccard – Age_binary)

Figure: PERMANOVA Jaccard for “Age_binary” (16S)

---

Fungi (ITS)

With Jaccard dissimilarity, establishment method (p = 0.001) and pH (p = 0.002) remained the strongest drivers, followed by site age (p = 0.006), ploughing (p = 0.024), and cutting (p = 0.04).

Pairwise comparisons showed the most significant differences between green hay/bush and both natural regeneration and seed mix sites (p = 0.001), and also between natural regeneration and seed mix (p = 0.004). Median Jaccard dissimilarity values were high across all establishment methods (0.75–0.85), indicating substantial variation in taxon presence/absence.

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough pH category Age category

View full QIIME2 results (Jaccard – Establishment)

Figure: PERMANOVA Jaccard for “Establishment” (ITS)

View full QIIME2 results (Jaccard – Cutting)

Figure: PERMANOVA Jaccard for “Cutting” (ITS)

View full QIIME2 results (Jaccard – Cattle)

Figure: PERMANOVA Jaccard for “Cattle” (ITS)

View full QIIME2 results (Jaccard – Sheep)

Figure: PERMANOVA Jaccard for “Sheep” (ITS)

View full QIIME2 results (Jaccard – Plough)

Figure: PERMANOVA Jaccard for “Plough” (ITS)

View full QIIME2 results (Jaccard – pH_category)

Figure: PERMANOVA Jaccard for “pH_category” (ITS)

View full QIIME2 results (Jaccard – Age_category)

Figure: PERMANOVA Jaccard for “Age_category” (ITS)

Unweighted Unifrac - Only for Bacteria (16S)

Emperor Plot

Here is a link to the Unweighted Unifrac Emperor Plot for more flexibility on QIIME2: Unweighted Unifrac Emperor Plot (16S)

Unweighted Unifrac Emperor Plot (16S)

PERMANOVA

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough Year_group pH binary Age binary

View full QIIME2 results (Unweighted Unifrac – Establishment)

Figure: PERMANOVA Unweighted Unifrac for “Establishment” (16S)

View full QIIME2 results (Unweighted Unifrac – Cutting)

Figure: PERMANOVA Unweighted Unifrac for “Cutting” (16S)

View full QIIME2 results (Unweighted Unifrac – Cattle)

Figure: PERMANOVA Unweighted Unifrac for “Cattle” (16S)

View full QIIME2 results (Unweighted Unifrac – Sheep)

Figure: PERMANOVA Unweighted Unifrac for “Sheep” (16S)

View full QIIME2 results (Unweighted Unifrac – Plough)

Figure: PERMANOVA Unweighted Unifrac for “Plough” (16S)

View full QIIME2 results (Unweighted Unifrac – Year_group)

Figure: PERMANOVA Unweighted Unifrac for “Year_group” (16S)

View full QIIME2 results (Unweighted Unifrac – pH_binary)

Figure: PERMANOVA Unweighted Unifrac for “pH_binary” (16S)

View full QIIME2 results (Unweighted Unifrac – Age_binary)

Figure: PERMANOVA Unweighted Unifrac for “Age_binary” (16S)

Weighted Unifrac - Only for Bacteria (16S)

Emperor Plot

Here is a link to the Weighted Unifrac Emperor Plot for more flexibility on QIIME2: Weighted Unifrac Emperor Plot (16S)

Weighted Unifrac Emperor Plot (16S)

PERMANOVA

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment Cutting Cattle Sheep Plough Year_group pH binary Age binary

View full QIIME2 results (Weighted Unifrac – Establishment)

Figure: PERMANOVA Weighted Unifrac for “Establishment” (16S)

View full QIIME2 results (Weighted Unifrac – Cutting)

Figure: PERMANOVA Weighted Unifrac for “Cutting” (16S)

View full QIIME2 results (Weighted Unifrac – Cattle)

Figure: PERMANOVA Weighted Unifrac for “Cattle” (16S)

View full QIIME2 results (Weighted Unifrac – Sheep)

Figure: PERMANOVA Weighted Unifrac for “Sheep” (16S)

View full QIIME2 results (Weighted Unifrac – Plough)

Figure: PERMANOVA Weighted Unifrac for “Plough” (16S)

View full QIIME2 results (Weighted Unifrac – Year_group)

Figure: PERMANOVA Weighted Unifrac for “Year_group” (16S)

View full QIIME2 results (Weighted Unifrac – pH_binary)

Figure: PERMANOVA Weighted Unifrac for “pH_binary” (16S)

View full QIIME2 results (Weighted Unifrac – Age_binary)

Figure: PERMANOVA Weighted Unifrac for “Age_binary” (16S)

---

Taxonomy Composition

Taxonomy Barplot

Bacteria (16S)

Microbial community composition varied notably across establishment methods :

• Green hay/bush soils were dominated by Alphaproteobacteria, Actinobacteria, and Bacilli

• Natural regeneration sites had higher abundances of Acidobacteria-6 and Thermoleophilia, suggesting adaptation to low-nutrient conditions

• Seed mix sites showed increased proportions of Alphaproteobacteria, Gammaproteobacteria, and Betaproteobacteria, indicating more copiotrophic (resource-demanding) communities typical of early-successional stages

Here is a link to the Taxonomy Barplots for more flexibility on QIIME2: Taxonomy Barplot (16S)

Select variable: Establishment Cutting Cattle Sheep Plough Age group pH binary Year group Age binary

Taxonomy Barplots Associated with 'Establishment' (16S)

Taxonomy Barplots Associated with 'Cutting' (16S)

Taxonomy Barplots Associated with 'Cattle' (16S)

Taxonomy Barplots Associated with 'Sheep' (16S)

Taxonomy Barplots Associated with 'Plough' (16S)

Taxonomy Barplots Associated with 'Age_group' (16S)

Taxonomy Barplots Associated with 'pH_binary' (16S)

Taxonomy Barplots Associated with 'Year_group' (16S)

Taxonomy Barplots Associated with 'Age_binary' (16S)

Fungi (ITS)

Here is a link to the Taxonomy Barplots for more flexibility on QIIME2: Taxonomy Barplot (ITS)

Select variable: Establishment Cutting Cattle Sheep Plough Age category pH category

Taxonomy Barplots Associated with 'Establishment' (ITS)

Taxonomy Barplots Associated with 'Cutting' (ITS)

Taxonomy Barplots Associated with 'Cattle' (ITS)

Taxonomy Barplots Associated with 'Sheep' (ITS)

Taxonomy Barplots Associated with 'Plough' (ITS)

Taxonomy Barplots Associated with 'Age_category' (ITS)

Taxonomy Barplots Associated with 'pH_category' (ITS)

Krona Plots

To explore the composition of soil microbial communities, we used Krona plots — interactive, circular charts that display taxonomic abundances in a hierarchical manner.

These plots allow users to intuitively navigate from broader taxonomic levels (such as Phylum) to more specific ones (like Genus), while simultaneously comparing relative abundances across taxa.

In this study, Krona plots provide a powerful and user-friendly way to:

• Visualise which microbial groups dominate each site
• Explore the taxonomic diversity present in bacterial, archaeal, and fungal communities

You can click on the images below to access the Krona plots for each site.

Bacteria (16S)

Krona Plot for Baltic_farm_1 (16S)

Fungi (ITS)

Krona Plot for Baltic_farm_1 (ITS)

Differential Abundance Analysis with ANCOM

We used ANCOM to identify taxa whose relative abundances significantly differed across groups. This method accounts for the compositional nature of microbiome data by comparing log-ratios between taxa. The results are shown as volcano-like plots, where the W statistic reflects how many pairwise comparisons a taxon was found to differ in. Significant taxa are highlighted accordingly.

Bacteria (16S)

ANCOM analysis revealed that while most bacterial taxa did not vary significantly across establishment methods, a subset—including Intrasporangiaceae and Blastococcus (both Actinobacteria)—showed strong differential abundance (W > 300, clr > 15). These taxa responded notably to restoration strategy, indicating that establishment method can drive specific shifts in bacterial community structure.

Similarly, the ANCOM volcano plot for pH showed that although most bacterial taxa remained stable across the pH gradient, several taxa were strongly pH-responsive (W > 200). These results suggest that soil pH exerts a selective influence, favouring distinct bacterial groups in either more acidic or alkaline conditions.

Select variable: Establishment Cutting Cattle Sheep Plough Age_group pH_binary Year_group Age_binary

View interactive volcano-like plot in QIIME2 (16S - Establishment )

Differences in Taxa Abundance Associated with 'Establishment' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Cutting )

Differences in Taxa Abundance Associated with 'Cutting' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Cattle )

Differences in Taxa Abundance Associated with 'Cattle' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Sheep )

Differences in Taxa Abundance Associated with 'Sheep' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Plough )

Differences in Taxa Abundance Associated with 'Plough' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Age_group )

Differences in Taxa Abundance Associated with 'Age_group' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - pH_binary )

Differences in Taxa Abundance Associated with 'pH_binary' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Year_group )

Differences in Taxa Abundance Associated with 'Year_group' (ANCOM Results) (16S)

View interactive volcano-like plot in QIIME2 (16S - Age_binary )

Differences in Taxa Abundance Associated with 'Age_binary' (ANCOM Results) (16S)

Fungi (ITS)

ANCOM analysis identified several fungal taxa as differentially abundant in relation to soil pH, site age, and establishment method:

• Metapochonia (W = 386) was significantly more abundant in acidic soils, suggesting a preference for low-pH environments.

• Two taxa were age-responsive: a Lasiosphaeriaceae taxon (W = 366) was more abundant in younger soils, while Paraphaeosphaeria (W = 320) was enriched in older soils, indicating shifts in fungal community structure with site maturation.

• Regarding establishment method, Gibellulopsis (W = 596) was predominantly found in green hay/bush and seed mix sites, but rare in natural regeneration sites. In contrast, Metarhizium (W = 535) was most abundant in naturally regenerating sites, likely reflecting adaptation to less-disturbed, self-organising ecosystems.

These findings highlight how restoration strategy, pH, and time since restoration can drive specific taxonomic responses, shaping fungal community composition in grassland soils.

Select variable: Establishment Cutting Cattle Sheep Plough Age_category pH_category

View interactive volcano-like plot in QIIME2 (ITS - Establishment )

Differences in Taxa Abundance Associated with 'Establishment' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - Cutting )

Differences in Taxa Abundance Associated with 'Cutting' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - Cattle )

Differences in Taxa Abundance Associated with 'Cattle' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - Sheep )

Differences in Taxa Abundance Associated with 'Sheep' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - Plough )

Differences in Taxa Abundance Associated with 'Plough' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - Age_category )

Differences in Taxa Abundance Associated with 'Age_category' (ANCOM Results) (ITS)

View interactive volcano-like plot in QIIME2 (ITS - pH_category )

Differences in Taxa Abundance Associated with 'pH_category' (ANCOM Results) (ITS)

---

Functional Diversity (Bacteria & Fungi)

Functional diversity refers to the variety of biological functions present in microbial communities. Unlike taxonomic analysis, which identifies who is there, functional analysis explores what they can do. This section includes two parallel analyses:

• Bacteria (16S): Functional profiles were inferred using PICRUSt2, predicting gene family and pathway abundance.
• Fungi (ITS): Functional roles were explored using FUNGuild, which assigns taxa to ecological guilds based on literature.

Alpha and beta diversity were analysed, along with differentially abundant functions or pathways.

Bacteria (16S)

Functional Alpha Diversity

The following metrics are based on PICRUSt2-inferred pathways.

We examined three complementary metrics: Shannon diversity (richness and evenness), Observed Features (richness only), and Evenness (dominance balance). You can use the menu below to explore how each metric varies according to different experimental variables.

Select variable: Establishment Cutting Cattle Sheep Plough Age_group pH_binary Age_binary

Shannon Index Boxplots

Here is a link to the Shannon Index Boxplot for more flexibility on QIIME2: Shannon Index Boxplot - Functional (16S)

Kruskal-Wallis p-value: 0.00164

Kruskal-Wallis p-value: 0.0496

Kruskal-Wallis p-value: 0.739

Kruskal-Wallis p-value: 0.681

Kruskal-Wallis p-value: 0.588

Kruskal-Wallis p-value: 0.144

Kruskal-Wallis p-value: 0.517

Kruskal-Wallis p-value: 0.0887

Observed Features Boxplots

Here is a link to the Observed features Boxplot for more flexibility on QIIME2: Observed Features Boxplot - Functional (16S)

Kruskal-Wallis p-value: 0.0354

Kruskal-Wallis p-value: 0.758

Kruskal-Wallis p-value: 0.969

Kruskal-Wallis p-value: 0.707

Kruskal-Wallis p-value: 0.923

Kruskal-Wallis p-value: 0.372

Kruskal-Wallis p-value: 0.292

Kruskal-Wallis p-value: 0.974

Evenness Significance Boxplots

Here is a link to the Evenness Boxplot for more flexibility on QIIME2: Evenness Boxplot - Functional (16S)

Kruskal-Wallis p-value: 0.0294

Kruskal-Wallis p-value: 0.13

Kruskal-Wallis p-value: 0.888

Kruskal-Wallis p-value: 0.858

Kruskal-Wallis p-value: 0.393

Kruskal-Wallis p-value: 0.84

Kruskal-Wallis p-value: 0.184

Kruskal-Wallis p-value: 0.235

Functional Beta Diversity

Emperor Plot – Bray-Curtis Distance

Based on predicted pathway abundance from PICRUSt2.

Here is a link to the Bray-Curtis Emperor Plot for more flexibility on QIIME2: Bray-Curtis Emperor Plot (16S)

Bray-Curtis Emperor Plot (Functional)

PERMANOVA – Bray-Curtis

Select variable: Establishment Cutting Cattle Sheep Plough Year group Age group Age Binary pH Binary

View full QIIME2 results (Bray–Curtis – Establishment)

Figure: PERMANOVA Bray-Curtis for “Establishment” (ITS)

View full QIIME2 results (Bray–Curtis – Cutting)

Figure: PERMANOVA Bray-Curtis for “Cutting” (ITS)

View full QIIME2 results (Bray–Curtis – Cattle)

Figure: PERMANOVA Bray-Curtis for “Cattle” (ITS)

View full QIIME2 results (Bray–Curtis – Sheep)

Figure: PERMANOVA Bray-Curtis for “Sheep” (ITS)

View full QIIME2 results (Bray–Curtis – Plough)

Figure: PERMANOVA Bray-Curtis for “Plough” (ITS)

View full QIIME2 results (Bray–Curtis – Year_group)

Figure: PERMANOVA Bray-Curtis for “Year_group” (ITS)

View full QIIME2 results (Bray–Curtis – Age_group)

Figure: PERMANOVA Bray-Curtis for “Age_group” (ITS)

View full QIIME2 results (Bray–Curtis – Age_binary)

Figure: PERMANOVA Bray-Curtis for “Age_binary” (ITS)

View full QIIME2 results (Bray–Curtis – pH_binary)

Figure: PERMANOVA Bray-Curtis for “pH_binary” (ITS)

Differentially Abundant Pathways

To identify pathways with significantly different abundances between groups, ANCOM was applied to PICRUSt2 pathway predictions.

Select variable: Establishment Cutting Cattle Sheep Plough Age_group pH_binary Age_binary

Differences in Pathway Abundance Associated with 'Establishment' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Cutting' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Cattle' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Sheep' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Plough' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Age_group' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'pH_binary' (ANCOM Results) (ITS)

Differences in Pathway Abundance Associated with 'Age_binary' (ANCOM Results) (ITS)

Heatmaps

These heatmaps display the relative abundance of predicted functional pathways across different samples, based on PICRUSt2-inferred metagenomic profiles from 16S rRNA gene data. Generated using QIIME2, each row represents a metabolic pathway (from MetaCyc), and each column corresponds to a sample group. The colour intensity indicates the predicted abundance of each functional pathway in the respective group.

You can click on each image to view a larger version in a new tab.

---

Fungi (ITS)

Functional Alpha Diversity

For fungal communities (ITS), functional diversity was assessed using metrics derived from predicted ecological functions (e.g., trophic modes, symbiosis potential).

The metrics below—Shannon diversity and Evenness—summarise the variability of ecological roles across different experimental conditions. Use the dropdown menu to explore patterns by variable.

Select variable: Establishment Cutting Cattle Sheep Plough Age category pH category Age

Shannon Index Boxplots

Here is a link to the Shannon Index Boxplot for more flexibility on QIIME2: Shannon Index Boxplot - Functional (ITS)

Kruskal-Wallis p-value: 0.529

Kruskal-Wallis p-value: 0.797

Kruskal-Wallis p-value: 0.133

Kruskal-Wallis p-value: 0.592

Kruskal-Wallis p-value: 0.161

Kruskal-Wallis p-value: 0.439

Kruskal-Wallis p-value: 0.00299

Kruskal-Wallis p-value: 0.372

Evenness Significance Boxplots

Here is a link to the Evenness Boxplot for more flexibility on QIIME2: Evenness Boxplot - Functional (ITS)

Kruskal-Wallis p-value: 0.226

Kruskal-Wallis p-value: 0.349

Kruskal-Wallis p-value: 0.106

Kruskal-Wallis p-value: 0.617

Kruskal-Wallis p-value: 0.331

Kruskal-Wallis p-value: 0.487

Kruskal-Wallis p-value: 0.00668

Kruskal-Wallis p-value: 0.206

Functional Beta Diversity

Bray-Curtis

Emperor Plot

Here is a link to the Bray-Curtis Emperor Plot for more flexibility on QIIME2: Bray-Curtis Emperor Plot (16S)

Bray-Curtis Emperor Plot (Fungi - ITS, Functional)

PERMANOVA

Select variable: Establishment pH category Age category

View full QIIME2 results (Bray–Curtis – Establishment)

Figure: PERMANOVA Bray-Curtis for “Establishment” (ITS)

View full QIIME2 results (Bray–Curtis – pH)

Figure: PERMANOVA Bray-Curtis for “pH” (ITS)

View full QIIME2 results (Bray–Curtis – Age)

Figure: PERMANOVA Bray-Curtis for “Age” (ITS)

Jaccard

Emperor Plot

Here is a link to the Jaccard Emperor Plot for more flexibility on QIIME2: Bray-Curtis Emperor Plot (16S)

Jaccard Emperor Plot (Fungi - ITS, Functional)

PERMANOVA

You can click on the images below to access the full QIIME2 report.

Select variable: Establishment pH category Age category

View full QIIME2 results (Jaccard – Establishment)

Figure: PERMANOVA Jaccard for “Establishment” (ITS)

View full QIIME2 results (Jaccard – pH)

Figure: PERMANOVA Jaccard for “pH” (ITS)

View full QIIME2 results (Jaccard – Age)

Figure: PERMANOVA Jaccard for “Age” (ITS)

Fungal Functional Guilds (FUNGuild)

In microbial ecology, a guild refers to a group of organisms that fulfil similar ecological roles, regardless of their taxonomic identity. Understanding functional guilds allows researchers to move beyond taxonomic profiles and assess the ecological roles that microbial communities may play in an environment.

To investigate the ecological roles of fungal communities, we used FUNGuild, a tool that assigns fungi to functional guilds based on curated databases and literature. These guilds represent ecological strategies such as:

• Saprotrophs: decomposers of organic matter
  
• Mycorrhizal fungi: symbionts associated with plant roots
  
• Pathogens: organisms that cause disease in plants or animals
• Ectomycorrhizal: symbionts that assist plant roots by enhancing nutrient and water uptake, especially nitrogen and phosphorus, while contributing to carbon cycling and soil nutrient mobilisation
• Arbuscular Mycorrhizal: symbionts that contribute to nutrient uptake
• Endophyte: microorganisms that promote the growth and development of the plants
• Lichenized: fungi that form symbiotic associations with algae or cyanobacteria, creating lichens that can survive in harsh environments by combining structural support and photosynthetic ability
• Parasites: fungi that live on or inside a host organism, extracting nutrients and often causing harm
• Symbiotrophs: fungi that engage in mutually beneficial relationships with host organisms, typically exchanging nutrients for resources like carbon or protection

This functional classification provides valuable insights into what fungi are likely doing in the ecosystem, beyond simply who they are.

This section explores the functional roles of fungi within each site, based on guild-level annotations provided by FUNGuild. Fungal guilds reflect ecological functions such as saprotrophy, symbiosis (e.g., mycorrhizal fungi), or pathogenicity. This approach provides insight into how fungal communities may contribute to ecosystem processes, complementing traditional taxonomic analyses.

Top 20 Most Abundant Guilds

The plot below highlights the top 20 most abundant fungal guilds identified using FUNGuild. To avoid clutter, the guild names are hidden on the y-axis; however, users can hover over each bar to reveal the full name, enabling interactive and detailed exploration of fungal functional diversity.

Top 20 functional guilds

Guild Abundance Across Sites

The figure below shows the total abundance of fungal ASVs across sites, aggregated by functional guild. This provides an overview of how guild-level composition varies between locations, which may reflect differences in land use, soil conditions, or restoration histories.

Fungal communities were dominated by Undefined Saprotrophs, Plant Pathogens, and Dung Saprotrophs. A high proportion of Unassigned fungi highlights limits in current taxonomic resolution.

Fungal Guild Abundance by Site

Guild Abundance by Variable

Establishment

Natural regeneration had the highest proportion of Undefined Saprotrophs, while Green Hay/Bush and Seed Mix had more Plant Pathogens and Dung Saprotrophs.

pH

Alkaline sites showed dominance of Undefined Saprotrophs; acidic sites had more Animal Parasites.

Age

Older sites had more Undefined Saprotrophs; younger sites showed higher proportions of Plant Pathogens and Dung Saprotrophs.

Select variable: Establishment Age pH

Fungal Guild Relative Abundance by pH_category

Fungal Guild Relative Abundance by Age_category

Fungal Guild Relative Abundance by Establishment

---

Correlation between Bacteria and Fungi

Taxonomic level

Functional level